#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#                                 Functions for validation                                ++
#                         Lei Zhu (leizhu@fas.harvard.edu), 11/08/2019                    ++
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Find common elments in two vectors
#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

common <- function(v1, v2){
  r1 <- rle(sort(v1))
  r2 <- rle(sort(v2))
  vals <- intersect(r1$values, r2$values)
  l1 <- r1$lengths[r1$values %in% vals]
  l2 <- r2$lengths[r2$values %in% vals]
  rep(vals, pmin(l1, l2))
}

#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Convert integer to character, with 0 heading
#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Integer2character	<- function(x){
  temp	<- as.character(x)
  if(x<10){
    temp	<- paste("0",temp,sep="")
  }
  return(temp)
}

#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Get seconds since 2000-01-01 00:00
#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

get_seconds_2010 <- function(c_date){
  second_temp <- as.integer(unclass(difftime(as.Date(c_date), as.Date("2010-01-01"), units = "days")))*24*3600
  return(second_temp)
}

#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Get AMF using GC and satellite profiles, with satellite scattering weighting functions
#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Cal_GC_AMF <- function( P_GC_pix , HCHO_GC_pix  , AD_GC_pix   , BXH_GC_pix, 
                        P_OMI_pix, Shape_OMI_pix, ScaW_OMI_pix, AMFg_OMI_pix ){
  
  # Get GC pressure at each level (center)
  P_TOP                             <- 1e-4 # Top of the atmosphere
  P_edge_GC_pix                     <- c(P_GC_pix,P_TOP)
  P_center_GC_pix                   <- array(NA,dim=c(NLEVS_GC))
  for(ilevel in 1:NLEVS_GC){
    P_center_GC_pix[ilevel]         <- (P_edge_GC_pix[ilevel] + P_edge_GC_pix[ilevel+1])*0.5
  }
  
  # Get GC HCHO number density
  HCHO_ND_GC_pix                    <- HCHO_GC_pix*1e-9*AD_GC_pix  # molec. cm-3
  
  # Project GC number density from GC grids to OMI grids
  GC_ND_OMI_grid                    <- splint(P_center_GC_pix, HCHO_ND_GC_pix, P_OMI_pix, df=length(HCHO_ND_GC_pix))
  
  # Get GC height vs P
  H_edge_GC_pix                     <- array(NA,dim=c(NLEVS_GC+1)) # km
  for(ilevel in 1:(NLEVS_GC+1)){
    if(ilevel==1){
      H_edge_GC_pix[ilevel]         <- 0
    }else{
      H_edge_GC_pix[ilevel]         <- H_edge_GC_pix[ilevel-1] + BXH_GC_pix[ilevel-1]*1e-3
    }
  }
  
  # Get OMI pressure at edge
  P_edge_OMI_pix                    <- array(NA,dim=c(NLEVS_GC+1))
  P_edge_OMI_pix[NLEVS_GC+1]        <- 0
  for(ilev in 1:NLEVS_GC){
    P_edge_OMI_pix[NLEVS_GC-ilev+1] <- P_OMI_pix[NLEVS_GC-ilev+1]*2 - P_edge_OMI_pix[NLEVS_GC-ilev+2]
  }
  
  # Get OMI height vs P
  H_edge_OMI                        <- splint(P_edge_GC_pix, H_edge_GC_pix, P_edge_OMI_pix, df=length(H_edge_GC_pix))
  # Can not higher than 80km
  H_edge_OMI[which(H_edge_OMI>80)]  <- 80
  
  # Compute sigma, dsigma, and height of each level
  OMI_edge_Sigma                    <- (P_edge_OMI_pix-min(P_edge_OMI_pix))/(max(P_edge_OMI_pix)-min(P_edge_OMI_pix))
  OMI_dSigma                        <- array(NA,dim=c(NLEVS_GC))
  OMI_dH                            <- array(NA,dim=c(NLEVS_GC))
  for(ilev in 1:NLEVS_GC){
    OMI_dSigma[ilev]                <- OMI_edge_Sigma[ilev] - OMI_edge_Sigma[ilev+1]
    OMI_dH[ilev]                    <- H_edge_OMI[ilev+1] - H_edge_OMI[ilev]
  }
  
  # Get shape factor
  GC_partial_column_OMI_grid        <- GC_ND_OMI_grid*OMI_dH*1e5 # molec. cm-2
  # Compute Obs Shape factor
  GC_ShapeFactor                    <- GC_partial_column_OMI_grid/sum(GC_partial_column_OMI_grid)/OMI_dSigma
  
  # Compute OMI Shape factor
  OMI_ShapeFactor                   <- Shape_OMI_pix/sum(Shape_OMI_pix)/OMI_dSigma
  
  # Compute AMF
  GC_AMF                            <- sum(ScaW_OMI_pix*GC_ShapeFactor*OMI_dSigma,na.rm=T)*AMFg_OMI_pix
  OMI_AMF                           <- sum(ScaW_OMI_pix*OMI_ShapeFactor*OMI_dSigma,na.rm=T)*AMFg_OMI_pix
  
  # Return
  return(c(GC_AMF, OMI_AMF))
}

#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Main validation script
#///////////////////////////////////////////////////////////////////////////////////////////
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Do_validation <- function(Campaign_ID){
  
  Campaign_NO             <- as.character(Campaign_info$No[Campaign_ID])
  Campaign_name           <- as.character(Campaign_info$Campaign[Campaign_ID])
  Campaign_date           <- seq(as.Date(Campaign_info$Date_Start[Campaign_ID]), as.Date(Campaign_info$Date_End[Campaign_ID]), "days")
  NDAYS                   <- length(Campaign_date)
  
  #---> Get plot region
  Plot_Lon_grid           <- seq(Campaign_info$PlotLon1[Campaign_ID]+0.5*Res,Campaign_info$PlotLon2[Campaign_ID]-0.5*Res,by=Res)
  Plot_Lat_grid           <- seq(Campaign_info$PlotLat1[Campaign_ID]+0.5*Res,Campaign_info$PlotLat2[Campaign_ID]-0.5*Res,by=Res)
  
  Plot_Domain_Lon         <- c(Campaign_info$PlotLon1[Campaign_ID],Campaign_info$PlotLon2[Campaign_ID],
                               Campaign_info$PlotLon2[Campaign_ID],Campaign_info$PlotLon1[Campaign_ID])
  Plot_Domain_Lat         <- c(Campaign_info$PlotLat2[Campaign_ID],Campaign_info$PlotLat2[Campaign_ID],
                               Campaign_info$PlotLat1[Campaign_ID],Campaign_info$PlotLat1[Campaign_ID])
  
  #---> Get study region
  Region_Lon              <- c(Campaign_info$StudyLon1[Campaign_ID],Campaign_info$StudyLon2[Campaign_ID],
                               Campaign_info$StudyLon2[Campaign_ID],Campaign_info$StudyLon1[Campaign_ID])
  Region_Lat              <- c(Campaign_info$StudyLat2[Campaign_ID],Campaign_info$StudyLat2[Campaign_ID],
                               Campaign_info$StudyLat1[Campaign_ID],Campaign_info$StudyLat1[Campaign_ID])
  
  #---> Get index
  Col_left                <- which.min(abs(min(Region_Lon)-Plot_Lon_grid))
  Col_right               <- which.min(abs(max(Region_Lon)-Plot_Lon_grid))
  Row_low                 <- which.min(abs(min(Region_Lat)-Plot_Lat_grid))
  Row_up                  <- which.min(abs(max(Region_Lat)-Plot_Lat_grid))
  
  print("")
  print(paste(rep("=",70),collapse = ""))
  print(paste("Doing validation for:", Campaign_NO, Campaign_name))
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                    Load flight data                                     ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  #---> Load flight RData
  RData_name	        <- paste(RData_folder,Campaign_NO,"_",Campaign_name,"_profile.RData",sep="")
  load(RData_name)
  
  #---> Load GC RData
  RData_name2	        <- paste(RData_folder,Campaign_NO,"_",Campaign_name,"_GC_column.RData",sep="")
  load(RData_name2)
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                      Plot 1                                             ++
  #                                   Flight track                                          ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if(Plot_Flight_tracks){
    
    dummy       <- array(-99,dim=c(length(Plot_Lon_grid),length(Plot_Lat_grid)))
    v_min       <- Flight_min
    v_max       <- Flight_max
    plot_height <- plot_height_all[Campaign_ID]
    plot_width  <- plot_width_all[Campaign_ID]
    
    plotname	<- paste(FigFolder,Campaign_NO,"_",Campaign_name,"_plot1_flight_tracks",".png",sep="")
    png(plotname, width=plot_width, height=plot_height, units="in", res=400)
    par(mar=c(2.5,1.7,1.5,1.5))# Bottom, left, top, right
    image.plot(Plot_Lon_grid,Plot_Lat_grid,dummy,zlim=c(v_min,v_max),
               horizontal=T,legend.shrink=1,axis.args = list(cex.axis =1,padj=-1.75,tck=0.2),
               legend.width=1,legend.mar=1,
               legend.args=list(text=expression(paste("HCHO [ppbv]")),padj=0.15,cex=1),          
               xlab='',ylab='',midpoint=T,axes=F,ann=F
    )
    title(xlab="",cex.lab=1.25,font.lab=2)
    axis(1,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,cex.axis=1.6,font=1,padj=-0.5)
    axis(3,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,labels=F)
    title(ylab="",cex.lab=1.25,font.lab=2)
    axis(2,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,cex.axis=1.6,font=1,las=1,hadj=0.5)
    axis(4,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,labels=F)
    title(main=paste(Campaign_NO,Campaign_name),cex.main=2,font.main=2)
    map('world',add=T,lwd=1.25,col="gray40")
    box(lwd=2)
    
    for(i_point in 1:Obs_points){
      col_id <- (Obs_MX[i_point,1]-v_min)/(v_max-v_min)*64
      if(col_id > 64){col_id <- 64}
      points(Obs_Lon[i_point],Obs_Lat[i_point],col=tim.colors()[col_id],cex=0.4,pch=16)
    }
    
    # Add polygon
    polygon(Region_Lon, Region_Lat, col = NULL, lty = 1, lwd = 2, border = "green")
    
    # Close plot
    dev.off()
    
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                      Plot 2                                             ++
  #                                   Vertical profiles                                     ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if(Plot_HCHO_profile){
    
    xmax 		  <- Profile_xmax
    ymax      <- Profile_ymax
    plot_y 		<- seq(0.5*(MAX_Height/Ngroups),(MAX_Height-0.5*(MAX_Height/Ngroups)),by=(MAX_Height/Ngroups))
    
    plotname	<- paste(FigFolder,Campaign_NO,"_",Campaign_name,"_plot2_profile_height.png",sep="")
    png(plotname, width=5, height=8, units="in", res=400)
    par(mar=c(4,4,3,3))# Bottom, left, top, right
    
    plot(Profile_GC_mean, plot_y, type="l",xlim=c(0, xmax),ylim=range(c(0,ymax)),axes=F,ann=F,lwd=1,col="red")
    axis(1, at=pretty(seq(0, xmax,by=0.5)),lwd=2,tck=0.02,padj=-0.8,cex.axis = 1.25)
    axis(2, at=seq(0,ymax,by=1),lwd=2,tck=0.02,hadj=0.4,cex.axis = 1.25,las=1)
    mtext(side = 2, text = "Altitude [km]" , line = 3, cex=2   , padj=0.75   )
    mtext(side = 1, text = "HCHO [ppb]"   , line = 3, cex=2   , padj=-0.75)
    mtext(side = 3, text = paste(Campaign_NO, Campaign_name), line = 3, cex=1.75, padj=1.5 )
    
    # add sd
    for(ii in 1:(length(Profile_Obs_mean[,1])-1)){
      y_temp <- c(plot_y[ii], plot_y[ii], plot_y[ii+1], plot_y[ii+1])
      x_temp <- c(Profile_Obs_mean[ii,1]-Profile_Obs_sd[ii],
                  Profile_Obs_mean[ii,1]+ Profile_Obs_sd[ii],
                  Profile_Obs_mean[ii+1,1]+ Profile_Obs_sd[ii+1],
                  Profile_Obs_mean[ii+1,1]-Profile_Obs_sd[ii+1])
      polygon(x_temp, y_temp,col="gray90",border="gray90")
    }
    
    lines(Profile_Obs_mean[,1], plot_y, type="l",xlim=c(0, xmax),ylim=range(c(0,MAX_Height)),lwd=3,col="black")
    lines(Profile_GC_mean, plot_y, type="l",xlim=c(0, xmax),ylim=range(c(0,MAX_Height)),lwd=3,col="red")
    
    text(2.25,6000,"Observation",col="black",cex=1.5)
    text(2.25,5500,"GEOS-Chem",col="red",cex=1.5)
    box(lwd=2)
    dev.off()
    
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                              Get GEOS-Chem and obsered columns                          ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  GC_Integrated_column   <- sum(Profile_GC_mean_P*100*6.02e23/8.314/Profile_GC_mean_T*1e-6*Profile_GC_mean*1e-9*500*1e2, na.rm=T)
  Obs_Integrated_column  <- sum(Profile_Obs_mean_P*100*6.02e23/8.314/Profile_Obs_mean_T*1e-6*Profile_Obs_mean[,1]*1e-9*500*1e2, na.rm=T)
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                 Read satellite L2 files                                 ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if(Read_HCHO_L2){
    

    #---> Get list of L2 files
    L2_file_list   <- list.files(L2_data_folder)
    L2_file_used   <- c()
 
    #---> Define data arraies
    Pixel_count    <- 0
    OMI_Date       <- array(NA,dim=c(MAX_pixels))
    OMI_Lat        <- array(NA,dim=c(MAX_pixels))
    OMI_Lon        <- array(NA,dim=c(MAX_pixels))
    OMI_VCD        <- array(NA,dim=c(MAX_pixels))
    OMI_SCD        <- array(NA,dim=c(MAX_pixels))
    OMI_Corr       <- array(NA,dim=c(MAX_pixels))
    OMI_VCD_new    <- array(NA,dim=c(MAX_pixels))
    OMI_AMF        <- array(NA,dim=c(MAX_pixels))
    GC_AMF         <- array(NA,dim=c(MAX_pixels))
    OMI_AMF2       <- array(NA,dim=c(MAX_pixels))
    OMI_AMFg       <- array(NA,dim=c(MAX_pixels))
    OMI_Albedo     <- array(NA,dim=c(MAX_pixels))
    OMI_GasProfile <- array(NA,dim=c(MAX_pixels,NLEVS_GC))
    OMI_ScatterW   <- array(NA,dim=c(MAX_pixels,NLEVS_GC))
    OMI_PreLevel   <- array(NA,dim=c(MAX_pixels,NLEVS_GC))
    
    #---> Read L2 files
    for(ifile in 1:length(L2_file_list)){
      
      date_L2      <- paste( substr(L2_file_list[ifile],20,23), substr(L2_file_list[ifile],25,26), substr(L2_file_list[ifile],27,28), sep="-")

      # Get day index
      DAY_temp     <- grep(date_L2, Campaign_date)
      
      # Is this date within the campaign period
      if( length(DAY_temp) > 0 ){
        
        L2_filename    <- L2_file_list[ifile]
        print(paste("  - Read", L2_filename))
        
        # Open hdf file and read
        file_name      <- paste(L2_data_folder,L2_filename,sep="")
        Time           <- h5read(file_name,Time_FN)
        Lat            <- h5read(file_name,Lat_FN)
        Lon            <- h5read(file_name,Lon_FN)
        SZA            <- h5read(file_name,SZA_FN)
        VZA            <- h5read(file_name,VZA_FN)
        CloudFra       <- h5read(file_name,CloudFra_FN)
        MainFlag       <- h5read(file_name,MainFlag_FN)
        GasProfile     <- h5read(file_name,GasProfile_FN) # molec. cm-2
        ScatterW       <- h5read(file_name,ScatterW_FN)   #
        PreLevel       <- h5read(file_name,PreLevel_FN)   # hPa
        VCD            <- h5read(file_name,VCD_FN)        #
        SCD            <- h5read(file_name,SCD_FN)        # later convert it to SCD
        AMF            <- h5read(file_name,AMF_FN)        #
        AMFg           <- h5read(file_name,AMFg_FN)       #
        Albedo         <- h5read(file_name,Albedo_FN)     #
        
        # Select the data
        NRows          <- dim(VCD)[1]
        NTracks        <- dim(VCD)[2]
        L2_file_used_1st <- 0
        
        # Any pixels in the domain?
        for(irow in 1:NRows){
          In_Domain    <- c()
          In_Domain    <- pnt.in.poly(cbind(Lon[irow,],Lat[irow,]),cbind(Plot_Domain_Lon, Plot_Domain_Lat))$pip
          if(sum(In_Domain)>0){
            L2_file_used_1st <- L2_file_used_1st + 1
            # Now need to save this file
            if(L2_file_used_1st==1){
              L2_file_used <- c(L2_file_used, L2_filename)
            }
            track_list <- which(In_Domain>0)
            for(itrack in 1:length(track_list)){
              if( MainFlag[irow, track_list[itrack]]  == 0            && 
                  VCD[irow, track_list[itrack]]       >= VCD_limit[1] && 
                  VCD[irow, track_list[itrack]]       <= VCD_limit[2] &&
                  CloudFra[irow, track_list[itrack]]  >= CF_limit[1]  && 
                  CloudFra[irow, track_list[itrack]]  <= CF_limit[2]  && 
                  SZA[irow, track_list[itrack]]       >= SZA_limit[1] && 
                  SZA[irow, track_list[itrack]]       <= SZA_limit[2] ){
                
                # Now get a valid pixel
                Pixel_count                  <- Pixel_count + 1
                OMI_Date[Pixel_count]        <- paste(Integer2character(Time[,track_list[itrack]][1]), Integer2character(Time[,track_list[itrack]][2]),
                                                      Integer2character(Time[,track_list[itrack]][3]),sep="")                    
                OMI_Lat[Pixel_count]         <- Lat[irow,track_list[itrack]]
                OMI_Lon[Pixel_count]         <- Lon[irow,track_list[itrack]]                            
                OMI_VCD[Pixel_count]         <- VCD[irow, track_list[itrack]]
                OMI_AMF[Pixel_count]         <- AMF[irow, track_list[itrack]]
                OMI_SCD[Pixel_count]         <- SCD[irow, track_list[itrack]]*AMF[irow, track_list[itrack]]
                OMI_Corr[Pixel_count]        <- SCD[irow, track_list[itrack]]*AMF[irow, track_list[itrack]] - 
                  VCD[irow, track_list[itrack]]*AMF[irow, track_list[itrack]]
                OMI_AMFg[Pixel_count]        <- AMFg[irow, track_list[itrack]]
                OMI_Albedo[Pixel_count]      <- Albedo[irow, track_list[itrack]]
                OMI_GasProfile[Pixel_count,] <- GasProfile[irow, track_list[itrack],]
                OMI_ScatterW[Pixel_count,]   <- ScatterW[irow, track_list[itrack],]/AMFg[irow, track_list[itrack]] # Need to consider AMFg
                OMI_PreLevel[Pixel_count,]   <- PreLevel[irow, track_list[itrack],]    
                
                # Re-compute the AMF using GC profile
                COL_GC_temp                  <- which.min(abs(OMI_Lon[Pixel_count]-Lon_grid_GC))[1]
                ROW_GC_temp                  <- which.min(abs(OMI_Lat[Pixel_count]-Lat_grid_GC))[1]
                P_GC_pix                     <- GC_P_1314LT_daily_avg[COL_GC_temp,ROW_GC_temp,,DAY_temp]
                HCHO_GC_pix                  <- GC_HCHO_1314LT_daily_avg[COL_GC_temp,ROW_GC_temp,,DAY_temp] # mixing ratio
                AD_GC_pix                    <- GC_AD_1314LT_daily_avg[COL_GC_temp,ROW_GC_temp,,DAY_temp]
                BXH_GC_pix                   <- GC_BXH_1314LT_daily_avg[COL_GC_temp,ROW_GC_temp,,DAY_temp]
                P_OMI_pix                    <- OMI_PreLevel[Pixel_count,]  
                Shape_OMI_pix                <- OMI_GasProfile[Pixel_count,]
                ScaW_OMI_pix                 <- OMI_ScatterW[Pixel_count,]
                AMFg_OMI_pix                 <- OMI_AMFg[Pixel_count]
                
                # Check valid GC output
                if(max(P_GC_pix)==0){
                  AMF_temp                   <- c(NA, NA)
                }else{
                  AMF_temp                   <- Cal_GC_AMF( P_GC_pix , HCHO_GC_pix  , AD_GC_pix   , BXH_GC_pix, 
                                                            P_OMI_pix, Shape_OMI_pix, ScaW_OMI_pix, AMFg_OMI_pix )  
                }
                
                GC_AMF[Pixel_count]          <- AMF_temp[1]
                # Compute AMF2 here to test 
                OMI_AMF2[Pixel_count]        <- AMF_temp[2]
                
                # Get updated VCD
                OMI_VCD_new[Pixel_count]     <- (OMI_SCD[Pixel_count]-OMI_Corr[Pixel_count])/GC_AMF[Pixel_count] 
              }	
            }	
          }
        } # Loop row (1-60)
      }
    }
    
    #---> Resize the arraies
    OMI_Date            <- OMI_Date[1:Pixel_count]
    OMI_Lat             <- OMI_Lat[1:Pixel_count]
    OMI_Lon             <- OMI_Lon[1:Pixel_count]
    OMI_VCD             <- OMI_VCD[1:Pixel_count]
    OMI_SCD             <- OMI_SCD[1:Pixel_count]
    OMI_Corr            <- OMI_Corr[1:Pixel_count]
    OMI_VCD_new         <- OMI_VCD_new[1:Pixel_count]
    OMI_AMF             <- OMI_AMF[1:Pixel_count]
    GC_AMF              <- GC_AMF[1:Pixel_count]
    OMI_AMF2            <- OMI_AMF2[1:Pixel_count]
    OMI_AMFg            <- OMI_AMFg[1:Pixel_count]
    OMI_Albedo          <- OMI_Albedo[1:Pixel_count]
    OMI_GasProfile      <- OMI_GasProfile[1:Pixel_count,]
    OMI_ScatterW        <- OMI_ScatterW[1:Pixel_count,]
    OMI_PreLevel        <- OMI_PreLevel[1:Pixel_count,]
    
    # Compare GC and OMI AMF
    reg                 <- lm(GC_AMF~OMI_AMF)
    summary(reg)
    
    # ------------------------------------------------------------------------------------------
    # Now regrid OMI pixels 
    # ------------------------------------------------------------------------------------------
    
    # For pixels within the dates with obs
    NCOLS_GRID          <- length(Plot_Lon_grid)
    NROWS_GRID          <- length(Plot_Lat_grid)
    Grid_N              <- array(0 ,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS ))
    
    # Save all pixels
    Grid_VCD_all        <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day))
    Grid_SCD_all        <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day))
    Grid_Corr_all       <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day))
    Grid_AMF_all        <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day))
    Grid_AMF_new_all    <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day))
    Grid_AMFg_all       <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day))
    Grid_Albedo_all     <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day))
    Grid_VCD_new_all    <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day))
    Grid_GasProfile_all <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day, NLEVS_GC))
    Grid_ScatterW_all   <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day, NLEVS_GC))
    Grid_PreLevel_all   <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, MAX_pixels_day, NLEVS_GC))
    
    # Daily Average
    Grid_VCD_d          <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS ))
    Grid_SCD_d          <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS))
    Grid_Corr_d         <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS))
    Grid_AMF_d          <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS))
    Grid_AMF_new_d      <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS))
    Grid_AMFg_d         <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS))
    Grid_Albedo_d       <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS))
    Grid_VCD_new_d      <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS ))
    Grid_GasProfile_d   <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, NLEVS_GC))
    Grid_ScatterW_d     <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, NLEVS_GC))
    Grid_PreLevel_d     <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NDAYS, NLEVS_GC))
    
    # Campigan averages
    Grid_VCD            <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID))
    Grid_SCD            <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID))
    Grid_Corr           <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID))
    Grid_AMF            <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID))
    Grid_AMF_new        <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID))
    Grid_AMFg           <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID))
    Grid_Albedo         <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID))
    Grid_VCD_new        <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID))
    Grid_GasProfile     <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NLEVS_GC))
    Grid_ScatterW       <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NLEVS_GC))
    Grid_PreLevel       <- array(NA,dim=c(NCOLS_GRID, NROWS_GRID, NLEVS_GC))
    
    #---> Loop pixels, store satellite values first
    for( ipixel in 1:Pixel_count ){
      Col_index         <- which.min(abs(Plot_Lon_grid-OMI_Lon[ipixel]))[1]
      Row_index         <- which.min(abs(Plot_Lat_grid-OMI_Lat[ipixel]))[1]
      
      # Is this pixel at the date with obs?
      Obs_flag          <- grep(paste(substr(OMI_Date[ipixel],1,4),substr(OMI_Date[ipixel],5,6),substr(OMI_Date[ipixel],7,8),sep="-"), Date_with_obs)
      if(length(Obs_flag)>0){
        Day_index        <- Date_with_obs_index[Obs_flag]
        Pixel_index                                                          <- Grid_N[Col_index, Row_index, Day_index] + 1
        Grid_N[Col_index, Row_index, Day_index]                              <- Pixel_index
        Grid_VCD_all[Col_index, Row_index, Day_index, Pixel_index ]          <- OMI_VCD[ipixel]	
        Grid_SCD_all[Col_index, Row_index, Day_index, Pixel_index ]          <- OMI_SCD[ipixel]
        Grid_Corr_all[Col_index, Row_index, Day_index, Pixel_index ]         <- OMI_Corr[ipixel]
        Grid_AMF_all[Col_index, Row_index, Day_index, Pixel_index ]          <- OMI_AMF[ipixel]	
        Grid_AMF_new_all[Col_index, Row_index, Day_index, Pixel_index ]      <- GC_AMF[ipixel]	
        Grid_AMFg_all[Col_index, Row_index, Day_index, Pixel_index ]         <- OMI_AMFg[ipixel]	
        Grid_Albedo_all[Col_index, Row_index, Day_index, Pixel_index ]       <- OMI_Albedo[ipixel]	
        Grid_VCD_new_all[Col_index, Row_index, Day_index, Pixel_index ]      <- OMI_VCD_new[ipixel]	
        Grid_GasProfile_all[Col_index, Row_index, Day_index, Pixel_index, ]  <- OMI_GasProfile[ipixel,]	
        Grid_ScatterW_all[Col_index, Row_index, Day_index, Pixel_index, ]    <- OMI_ScatterW[ipixel,]	
        Grid_PreLevel_all[Col_index, Row_index, Day_index, Pixel_index, ]    <- OMI_PreLevel[ipixel,]		
      }
    }
    
    #---> Do regridding
    for(icol in 1: NCOLS_GRID){
      for(irow in 1: NROWS_GRID){
        
        # Get daily averages first
        for(iday in 1: NDAYS){
          if(Grid_N[icol, irow, iday] > Pixel_limit){
            Grid_VCD_d[icol, irow, iday]      <- mean(Grid_VCD_all[icol, irow, iday, ],na.rm=T)
            Grid_SCD_d[icol, irow, iday]      <- mean(Grid_SCD_all[icol, irow, iday, ],na.rm=T)
            Grid_Corr_d[icol, irow, iday]     <- mean(Grid_Corr_all[icol, irow, iday, ],na.rm=T)
            Grid_AMF_d[icol, irow, iday]      <- mean(Grid_AMF_all[icol, irow, iday, ],na.rm=T)
            Grid_AMF_new_d[icol, irow, iday]  <- mean(Grid_AMF_new_all[icol, irow, iday, ],na.rm=T)
            Grid_AMFg_d[icol, irow, iday]     <- mean(Grid_AMFg_all[icol, irow, iday, ],na.rm=T)			
            Grid_Albedo_d[icol, irow, iday]   <- mean(Grid_Albedo_all[icol, irow, iday, ],na.rm=T)			
            Grid_VCD_new_d[icol, irow, iday]  <- mean(Grid_VCD_new_all[icol, irow, iday, ],na.rm=T)
            for(ilayer in 1:NLEVS_GC){
              Grid_GasProfile_d[icol, irow, iday, ilayer] <- mean(Grid_GasProfile_all[icol, irow, iday, , ilayer],na.rm=T)
              Grid_ScatterW_d[icol, irow, iday, ilayer]   <- mean(Grid_ScatterW_all[icol, irow, iday, , ilayer],na.rm=T)
              Grid_PreLevel_d[icol, irow, iday, ilayer]   <- mean(Grid_PreLevel_all[icol, irow, iday, , ilayer],na.rm=T)				
            }
          }
        }
        
        # Now average all obs based on daily averages
        Grid_VCD[icol, irow]      <- mean(Grid_VCD_d[icol, irow, ],na.rm=T)
        Grid_SCD[icol, irow]      <- mean(Grid_SCD_d[icol, irow, ],na.rm=T)
        Grid_Corr[icol, irow]     <- mean(Grid_Corr_d[icol, irow,],na.rm=T)
        Grid_AMF[icol, irow]      <- mean(Grid_AMF_d[icol, irow, ],na.rm=T)
        Grid_AMF_new[icol, irow]  <- mean(Grid_AMF_new_d[icol, irow, ],na.rm=T)
        Grid_AMFg[icol, irow]     <- mean(Grid_AMFg_d[icol, irow, ],na.rm=T)			
        Grid_Albedo[icol, irow]   <- mean(Grid_Albedo_d[icol, irow, ],na.rm=T)
        Grid_VCD_new[icol, irow]  <- mean(Grid_VCD_new_d[icol, irow, ],na.rm=T)
        for(ilayer in 1:NLEVS_GC){
          Grid_GasProfile[icol, irow, ilayer] <- mean(Grid_GasProfile_d[icol, irow, , ilayer],na.rm=T)
          Grid_ScatterW[icol, irow, ilayer]   <- mean(Grid_ScatterW_d[icol, irow, , ilayer],na.rm=T)
          Grid_PreLevel[icol, irow, ilayer]   <- mean(Grid_PreLevel_d[icol, irow, , ilayer],na.rm=T)				
        }
        
      }
    }
    
    # ------------------------------------------------------------------------------------------
    # Save RData
    # ------------------------------------------------------------------------------------------
    RData_name5	          <- paste(RData_folder,Campaign_NO,"_",Campaign_name,"_satellite_column.RData",sep="")
    
    save(file=RData_name5, Grid_GasProfile, Grid_ScatterW     , Grid_PreLevel, Plot_Lat_grid, Plot_Lon_grid,
         Grid_VCD_d , Grid_VCD_new_d, Grid_VCD        , Grid_SCD , Grid_VCD_new , Grid_Corr    , 
         Grid_AMF   , Grid_AMFg         , Grid_Albedo, L2_file_used )
    
    write.csv2(L2_file_used, row.names=F, file=paste(RData_folder,Campaign_NO,"_",Campaign_name,"_satellite_list",sep=""))
    
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                  Load satellite data                                    ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  RData_name5	        <- paste(RData_folder,Campaign_NO,"_",Campaign_name,"_satellite_column.RData",sep="")
  load(RData_name5)
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                       Get mean profiles and scattering weights                          ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  #---> Get mean a-prior profile, pressure levels, and scattering weights
  #     Averaged over the domain during the campigan
  OMI_GasProfile_avg             <- array(NA,dim=c(47))
  OMI_ScatterW_avg               <- array(NA,dim=c(47))
  OMI_PreLevel_avg               <- array(NA,dim=c(47))
  OMI_AMFg_avg                   <- c()
  
  for(ilayer in 1:47){
    OMI_GasProfile_avg[ilayer] <- mean(Grid_GasProfile[Col_left:Col_right, Row_low:Row_up, ilayer], na.rm=T)
    OMI_ScatterW_avg[ilayer]   <- mean(Grid_ScatterW[Col_left:Col_right  , Row_low:Row_up, ilayer], na.rm=T)
    OMI_PreLevel_avg[ilayer]   <- mean(Grid_PreLevel[Col_left:Col_right  , Row_low:Row_up, ilayer], na.rm=T)	
  }
  
  OMI_AMFg_avg                   <- mean(Grid_AMFg[Col_left:Col_right  , Row_low:Row_up], na.rm=T)
  
  # ------------------------------------------------------------------------------------------
  # Convert Pressure to Sigma coordinates
  # ------------------------------------------------------------------------------------------
  
  #---> Compute the pressure at edge
  OMI_PreLevel_edge <- array(NA,dim=c(48))
  OMI_PreLevel_edge[48]          <- 0
  for(ilayer in 1:47){
    OMI_PreLevel_edge[48-ilayer] <- OMI_PreLevel_avg[48-ilayer]*2 - OMI_PreLevel_edge[48-ilayer+1]
  }
  
  #---> Compute sigma, dsigma, and height of each level
  OMI_Sigma_center               <- array(NA,dim=c(47))
  OMI_dSigma                     <- array(NA,dim=c(47))
  OMI_dH                         <- array(NA,dim=c(47))
  
  OMI_Sigma_edge                 <- (OMI_PreLevel_edge-min(OMI_PreLevel_edge))/(max(OMI_PreLevel_edge)-min(OMI_PreLevel_edge))
  
  #---> Get the height of each level
  OMI_H_edge                     <- splint(Standard_Air_P, Standard_Air_H, OMI_PreLevel_edge, df=length(Standard_Air_H))
  
  #---> Can not higher than 80km
  OMI_H_edge[which(OMI_H_edge>80)] <- 80
  for(ilayer in 1:47){
    OMI_dSigma[ilayer]           <- OMI_Sigma_edge[ilayer] - OMI_Sigma_edge[ilayer+1]
    OMI_Sigma_center[ilayer]     <- (OMI_Sigma_edge[ilayer] + OMI_Sigma_edge[ilayer+1])*0.5
    OMI_dH[ilayer]               <- OMI_H_edge[ilayer+1] - OMI_H_edge[ilayer]
  }
  
  #---> Compute OMI Shape factor
  OMI_ShapeFactor_avg              <- OMI_GasProfile_avg/sum(OMI_GasProfile_avg)/OMI_dSigma
  
  #---> Load ATom-1 and 2 obs profiles, use for missing observations above
  if( Campaign_ID!=11 && Campaign_ID!=12){
    Profile_Obs_mean_ATom            <- (Profile_Obs_mean_ATom1 + Profile_Obs_mean_ATom2)/2 # ppb
    Profile_Obs_mean_P_ATom          <- (Profile_Obs_mean_P_ATom1 + Profile_Obs_mean_P_ATom2)/2
    Profile_Obs_mean_T_ATom          <- (Profile_Obs_mean_T_ATom1 + Profile_Obs_mean_T_ATom2)/2
    
    for(ilevel in 1:length(Profile_Obs_mean[,1])){
      if( is.na(Profile_Obs_mean[ilevel,1]) && !is.na(Profile_Obs_mean_ATom[ilevel]) ){
        Profile_Obs_mean[ilevel,1] <- Profile_Obs_mean_ATom[ilevel]
        Profile_Obs_mean_P[ilevel] <- Profile_Obs_mean_P_ATom[ilevel]
        Profile_Obs_mean_T[ilevel] <- Profile_Obs_mean_T_ATom[ilevel]
      }
    }
  }
  
  Obs_number_desnity               <- Profile_Obs_mean_P*100*6.02e23/8.314/Profile_Obs_mean_T*1e-6*Profile_Obs_mean[,1]*1e-9 # molec. cm-3
  
  #---> Make sure no extrapolation
  Obs_number_desnity_OMI_grid      <- array(NA,dim=c(47))
  OMI_lev1                         <- min(which( (OMI_PreLevel_avg-max(Profile_Obs_mean_P, na.rm=T))<0 ))
  OMI_lev2                         <- max(which( (OMI_PreLevel_avg-min(Profile_Obs_mean_P, na.rm=T))>0 ))
  Profile_Obs_mean_P2              <- Profile_Obs_mean_P[which(!is.na(Profile_Obs_mean_P))]
  Obs_number_desnity2              <- Obs_number_desnity[which(!is.na(Obs_number_desnity))]
  
  #---> Project number density and T to the OMI grids
  Obs_number_desnity_OMI_grid[OMI_lev1:OMI_lev2] <- splint(Profile_Obs_mean_P2, Obs_number_desnity2, OMI_PreLevel_avg[OMI_lev1:OMI_lev2], df=length(Obs_number_desnity2))
  Obs_partial_density_OMI_grid     <- Obs_number_desnity_OMI_grid*OMI_dH*1e5 #m olec. cm-2
  
  #---> Fill the NA with OMI values
  Obs_partial_density_OMI_grid[which(is.na(Obs_partial_density_OMI_grid))] <- OMI_GasProfile_avg[which(is.na(Obs_partial_density_OMI_grid))]
  
  #---> Compute Obs Shape factor
  Obs_ShapeFactor_avg              <- Obs_partial_density_OMI_grid/sum(Obs_partial_density_OMI_grid)/OMI_dSigma
  
  #---> Compute AMF
  AMF_OMI                          <- sum(OMI_ScatterW_avg*OMI_ShapeFactor_avg*OMI_dSigma,na.rm=T)*mean(OMI_AMFg_avg,na.rm=T)
  AMF_Obs                          <- sum(OMI_ScatterW_avg*Obs_ShapeFactor_avg*OMI_dSigma,na.rm=T)*mean(OMI_AMFg_avg,na.rm=T)
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                      Plot 3                                             ++
  #                         OMI profiles and scattering weights                             ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if(Plot_AMF_components){
    
    plotname	<- paste(FigFolder,Campaign_NO,"_",Campaign_name,"_plot3_AMF.png",sep="")
    
    png(plotname, width=6, height=6, units="in", res=400)
    par(mar = c(4.5, 5, 5.5, 5) ) # Leave space for z axis
    x_min <- 0
    x_max <- 10
    plot_pressure <- rev(seq(0,1000,by=50))
    plot(OMI_ShapeFactor_avg, OMI_PreLevel_avg,type="l",axes=F,bty="n",xlab="",ylab="",col="blue",xlim=c(x_min,x_max),lwd=2,lty=1,ylim=rev(range(OMI_PreLevel_avg)))
    lines(Obs_ShapeFactor_avg, OMI_PreLevel_avg,type="l",bty="n",xlab="",ylab="",col="black",xlim=c(x_min,x_max),lwd=2,lty=1,ylim=rev(range(OMI_PreLevel_avg)))
    
    text(0.6*x_max,150,paste("Observed (AMF=",format(round(AMF_Obs, 2), nsmall = 2), ")",sep=""),col="black",cex=1.25)
    text(0.6*x_max,220,paste("OMI (AMF="     ,format(round(AMF_OMI, 2), nsmall = 2), ")",sep=""),col="blue" ,cex=1.25)
    
    axis(side=1, at = pretty(seq(x_min,x_max,by=0.1)),cex.axis=1.25, lwd=2, padj=-0.7, tck=-0.01)
    
    par(new = TRUE)
    x_min2 <- 0
    x_max2 <- 1.2
    plot(OMI_ScatterW_avg, OMI_PreLevel_avg,type="l",axes=FALSE,bty="n",xlab="",ylab="",col="blue",xlim=c(x_min2,x_max2),lwd=2,lty=5,ylim=rev(range(OMI_PreLevel_avg))) 
    lines(OMI_ScatterW_avg*OMI_ShapeFactor_avg, OMI_PreLevel_avg,type="l",bty="n",xlab="",ylab="",col="blue",xlim=c(x_min2,x_max2),lwd=2,lty=3,ylim=rev(range(OMI_PreLevel_avg))) 
    lines(OMI_ScatterW_avg*Obs_ShapeFactor_avg, OMI_PreLevel_avg,type="l",bty="n",xlab="",ylab="",col="black",xlim=c(x_min2,x_max2),lwd=2,lty=3,ylim=rev(range(OMI_PreLevel_avg))) 
    
    axis(side=2, at = pretty(rev(c(0,1000))),cex.axis=1.25, las=2, hadj=0.7, lwd=2, tck=-0.01)
    axis(side=3, at = c(0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4),cex.axis=1.25, padj=0.7, lwd=2, tck=-0.01)
    
    mtext(expression(paste("Pressure [hPa]")), side=2, padj=-0.5, line=2,cex=1.5)
    mtext(expression(paste("S(",sigma,") (solid)")), side=1, line=2, padj=0.2, cex=1.5)
    mtext(expression(paste("w(",sigma,") (dashed); ", "S(",sigma,")w(",sigma,") (dotted)" )), side=3, line=2,padj=0.4,cex=1.5)
    
    par(new = TRUE)
    plot(-99,-99,type="l",axes=FALSE,bty="n",xlab="",ylab="",col="black",xlim=c(x_min2,x_max2),lwd=4,lty=5,ylim=rev(range(seq(0,1,by=0.1))))  
    axis(side=4, pretty(rev(range(seq(0,1,by=0.1)))),cex.axis=1.25, las=2, hadj=0.5, lwd=2, tck=-0.01)
    mtext(expression(paste("Sigma coordinate (", sigma, ")")), side=4, line=2, padj=0.4, cex=1.5)
    
    box(lwd=3)
    dev.off()
    
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                      Plot 4                                             ++
  #                               Satellite mean columns                                    ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if(Plot_Satellite_column){
    
    zmin        <- Satellite_min
    zmax        <- Satellite_max[Campaign_ID]
    
    plotname	<- paste(FigFolder,Campaign_NO,"_",Campaign_name,"_plot4_OMI_column",".png",sep="")
    png(plotname, width=plot_width, height=plot_height, units="in", res=400)
    par(mar=c(2.5,1.5,1.5,1))# Bottom, left, top, right
    image.plot(Plot_Lon_grid,Plot_Lat_grid, Grid_VCD/1e15,zlim=c(zmin,zmax),
               horizontal=T,legend.shrink=1,axis.args = list(cex.axis =1,padj=-1.75,tck=0.2),
               legend.width=1,legend.mar=1,
               legend.args=list(text=expression(paste("HCHO column [",10^15," molecules ",cm^-2,"]",sep="")),padj=0.15,cex=1),           
               xlab='',ylab='',midpoint=T,axes=F,ann=F
    )
    title(xlab="",cex.lab=1.25,font.lab=2)
    axis(1,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,cex.axis=1,font=1,padj=-1.5)
    axis(3,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,labels=F)
    title(ylab="",cex.lab=1.25,font.lab=2)
    axis(2,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,cex.axis=1,font=1,las=1,hadj=0.2)
    axis(4,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,labels=F)
    
    title(main=paste(Campaign_NO,Campaign_name,"OMI column"),cex.main=1.1,font.main=2)
    map('world',add=T,lwd=1.25,col="gray40")
    polygon(Region_Lon, Region_Lat, col = NULL, lty = 1, lwd = 3, border = "green")
    box(lwd=2)
    dev.off()
    
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                      Plot 5                                             ++
  #                        GEOS-Chem mean column, corrected                                 ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if(Plot_Corr_GC_column){
    
    zmin        <- CorrectedGC_min
    zmax        <- CorrectedGC_max[Campaign_ID]
    
    plotname            <- paste(FigFolder,Campaign_NO,"_",Campaign_name,"_plot5_Corrected_GC_column",".png",sep="")
    png(plotname, width=plot_width, height=plot_height, units="in", res=400)
    par(mar=c(2.5,1.5,1.5,1))# Bottom, left, top, right
    image.plot(Plot_Lon_grid,Plot_Lat_grid, GC_HCHO_OMI_grids$z/1e15*Obs_Integrated_column/GC_Integrated_column,zlim=c(zmin, zmax),
               horizontal=T,legend.shrink=1,axis.args = list(cex.axis =1,padj=-1.75,tck=0.2),
               legend.width=1,legend.mar=1,
               legend.args=list(text=expression(paste("HCHO column [",10^15," molecules ",cm^-2,"]",sep="")),padj=0.15,cex=1),           
               xlab='',ylab='',midpoint=T,axes=F,ann=F
    )
    title(xlab="",cex.lab=1.25,font.lab=2)
    axis(1,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,cex.axis=1,font=1,padj=-1.5)
    axis(3,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,labels=F)
    title(ylab="",cex.lab=1.25,font.lab=2)
    axis(2,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,cex.axis=1,font=1,las=1,hadj=0.2)
    axis(4,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,labels=F)
    
    title(main=paste(Campaign_NO,Campaign_name,"corrected GC column"),cex.main=1.1,font.main=2)
    map('world',add=T,lwd=1.25,col="gray40")
    polygon(Region_Lon, Region_Lat, col = NULL, lty = 1, lwd = 3, border = "green")
    box(lwd=2)
    dev.off()
    
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                      Plot 6                                             ++
  #                           Updated Satellite mean columns                                ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if(Plot_Updated_column){
    
    zmin        <- NewSatellite_min
    zmax        <- NewSatellite_max[Campaign_ID]
    
    plotname	<- paste(FigFolder,Campaign_NO,"_",Campaign_name,"_plot6_Updated_column",".png",sep="")
    png(plotname, width=plot_width, height=plot_height, units="in", res=400)
    par(mar=c(2.5,1.5,1.5,1))# Bottom, left, top, right
    image.plot(Plot_Lon_grid,Plot_Lat_grid, Grid_VCD_new/1e15,zlim=c(zmin,zmax),
               horizontal=T,legend.shrink=1,axis.args = list(cex.axis =1,padj=-1.75,tck=0.2),
               legend.width=1,legend.mar=1,
               legend.args=list(text=expression(paste("HCHO column [",10^15," molecules ",cm^-2,"]",sep="")),padj=0.15,cex=1),           
               xlab='',ylab='',midpoint=T,axes=F,ann=F
    )
    title(xlab="",cex.lab=1.25,font.lab=2)
    axis(1,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,cex.axis=1,font=1,padj=-1.5)
    axis(3,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,labels=F)
    title(ylab="",cex.lab=1.25,font.lab=2)
    axis(2,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,cex.axis=1,font=1,las=1,hadj=0.2)
    axis(4,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,labels=F)
    
    title(main=paste(Campaign_NO,Campaign_name,"updated column"),cex.main=1.1,font.main=2)
    map('world',add=T,lwd=1.25,col="gray40")
    polygon(Region_Lon, Region_Lat, col = NULL, lty = 1, lwd = 3, border = "green")
    box(lwd=2)
    dev.off()
    
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                                      Plot 7                                             ++
  #                                      Albedo                                             ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  if(Plot_Albedo){
    
    zmin        <- Albedo_min
    zmax        <- Albedo_max
    
    plotname	<- paste(FigFolder,Campaign_NO,"_",Campaign_name,"_plot7_Albedo",".png",sep="")
    png(plotname, width=plot_width, height=plot_height, units="in", res=400)
    par(mar=c(2.5,1.5,1.5,1))# Bottom, left, top, right
    image.plot(Plot_Lon_grid,Plot_Lat_grid, Grid_Albedo,zlim=c(zmin,zmax),
               horizontal=T,legend.shrink=1,axis.args = list(cex.axis =1,padj=-1.75,tck=0.2),
               legend.width=1,legend.mar=1,
               legend.args=list(text=expression(paste("Surface Albedo [unitless]",sep="")),padj=0.15,cex=1),           
               xlab='',ylab='',midpoint=T,axes=F,ann=F
    )
    title(xlab="",cex.lab=1.25,font.lab=2)
    axis(1,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,cex.axis=1,font=1,padj=-1.5)
    axis(3,at=pretty(Plot_Lon_grid),tck=0.015,lwd=1.5,labels=F)
    title(ylab="",cex.lab=1.25,font.lab=2)
    axis(2,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,cex.axis=1,font=1,las=1,hadj=0.2)
    axis(4,at=pretty(Plot_Lat_grid),tck=0.015,lwd=1.5,labels=F)
    
    title(main=paste(Campaign_NO,Campaign_name,"Surface Albedo"),cex.main=1.1,font.main=2)
    map('state',add=T,lwd=1.25,col="gray40")
    polygon(Region_Lon, Region_Lat, col = NULL, lty = 1, lwd = 3, border = "green")
    box(lwd=2)
    dev.off()
    
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                           Get the mean values within the doamin                         ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  Satellite_mean_VCD            <- mean(Grid_VCD[Col_left: Col_right, Row_low: Row_up],na.rm=T)
  Satellite_mean_AMF            <- mean(Grid_AMF[Col_left: Col_right, Row_low: Row_up],na.rm=T)
  Satellite_mean_SCD            <- mean(Grid_SCD[Col_left: Col_right, Row_low: Row_up],na.rm=T)
  Satellite_mean_Corr           <- mean(Grid_Corr[Col_left: Col_right, Row_low: Row_up],na.rm=T)
  Satellite_mean_VCD1           <- mean((Grid_SCD[Col_left: Col_right, Row_low: Row_up]-
                                           Grid_Corr[Col_left: Col_right, Row_low: Row_up])/
                                          Grid_AMF[Col_left: Col_right, Row_low: Row_up],na.rm=T)
  Satellite_mean_VCD2           <- (Satellite_mean_SCD-Satellite_mean_Corr)/Satellite_mean_AMF
  
  Updated_satellite_mean_column <- (Satellite_mean_SCD-Satellite_mean_Corr)/AMF_Obs
  GC_mean_column                <- mean(GC_HCHO_OMI_grids$z[Col_left: Col_right, Row_low: Row_up],na.rm=T)
  Obs_informed_GC_mean_column   <- GC_mean_column/( GC_Integrated_column/Obs_Integrated_column)
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                    Get spatial corraltion between columns and albedo                    ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  Satellite_SCD     <- c() 
  Satellite_Albedo1 <- c()
  
  for(col in Col_left: Col_right){
    for(row in Row_low: Row_up){
      if( !is.na(Grid_SCD[col,row]) && !is.na(Grid_Albedo[col,row]) ){
        Satellite_SCD     <- c(Satellite_SCD, Grid_SCD[col,row])
        Satellite_Albedo1 <- c(Satellite_Albedo1,Grid_Albedo[col,row] )
      }
    }
  }
  
  Satellite_VCD    <- c() 
  Satellite_Albedo2 <- c()
  for(col in Col_left: Col_right){
    for(row in Row_low: Row_up){
      if( !is.na(Grid_VCD[col,row]) && !is.na(Grid_Albedo[col,row]) ){
        Satellite_VCD     <- c(Satellite_VCD, Grid_VCD[col,row])
        Satellite_Albedo2 <- c(Satellite_Albedo2,Grid_Albedo[col,row] )
      }
    }
  }
  
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #                              Write validation results to the file                       ++
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  cat("", file=logFile, append=TRUE, sep = "\n")
  cat(paste(rep("=",60),collapse = ""), file=logFile, append=TRUE, sep = "\n")
  cat(paste(Campaign_NO, Campaign_name), file=logFile, append=TRUE, sep = "\n")
  cat(paste(rep("=",60),collapse = ""), file=logFile, append=TRUE, sep = "\n")
  
  cat(paste(" ==> GC mean column (1E+15)               : ", round(GC_mean_column/1e15,2) ), file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Obs-informed GC mean column (1E+15)  : ", round(Obs_informed_GC_mean_column/1e15,2) ), file=logFile, append=TRUE, sep = "\n")
  
  cat("", file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Mean satellite AMFg                  : ", round(OMI_AMFg_avg,2 )), file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Mean satellite AMF                   : ", round(Satellite_mean_AMF,2) ), file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Satellite AMF using mean profiles    : ", round(AMF_OMI,2) ), file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Mean SCD                             : ", round(Satellite_mean_SCD/1e15,2) ), file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Mean correction                      : ", round(Satellite_mean_Corr/1e15,2) ), file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Satellite mean VCD, averaged (1E+15) : ", round(Satellite_mean_VCD/1e15,2)  ), file=logFile, append=TRUE, sep = "\n")
  cat(paste("     Mean bias (%)                        : ", format(round(Satellite_mean_VCD/Obs_informed_GC_mean_column-1,3)*100,0)), file=logFile, append=TRUE, sep = "\n")
  
  cat("", file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Satellite mean VCD, computed (1E+15) : ", round(Satellite_mean_VCD2/1e15,2)  ), file=logFile, append=TRUE, sep = "\n")
  cat(paste("     Mean bias (%)                        : ", format(round(Satellite_mean_VCD2/Obs_informed_GC_mean_column-1,3)*100,0)), file=logFile, append=TRUE, sep = "\n")
  
  cat("", file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Mean AMF using observed profiles     : ", round(AMF_Obs,2) ), file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Updated satellite column (1E+15)     : ", round(Updated_satellite_mean_column/1e15,2) ), file=logFile, append=TRUE, sep = "\n")
  cat(paste("     Mean bias (%)                        : ", format(round(Updated_satellite_mean_column/Obs_informed_GC_mean_column-1,3)*100,0)), file=logFile, append=TRUE, sep = "\n")
  
  cat("", file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Correlation between SCD and Albedo   : ", round(cor(Satellite_SCD,Satellite_Albedo1),3 )), file=logFile, append=TRUE, sep = "\n")
  cat(paste(" ==> Correlation between VCD and Albedo   : ", round(cor(Satellite_VCD,Satellite_Albedo2),3) ), file=logFile, append=TRUE, sep = "\n")

}
